# Dickstein, Ho, and Mark (2023)
# This script collects results for the bootstraped counterfactuals

# Preliminaries
setwd("../../../library")
source("PreliminariesCode.R")
setwd(paste0(project_folder, "/analysis/counterfactuals/bootstrap"))

## Which spec do you want to use?
spectouse <- "lowrisk_simplemh_censor0.2"

## Where is the counterfactual folder?
counterfactual_folder <- paste0(project_folder, "/analysis/counterfactuals")

# What alpha-psi cutoff to use: 
alphampsicutoff <- .05

# What is the vector of bootstrap draws?
k_vec <- 1:20

# Where is the output folder
output_folder <-  paste0(counterfactual_folder, "/bootstrap/output")

# Load data

# Counterfactual Results
results_list_base <- list()
results_list_1 <- list()
results_list_3 <- list()
results_list_5 <- list()
for(i in k_vec){
 print(paste0("LOADING COUNTERFACTUAL", i))
  load(paste0(output_folder, "/output_base_results_draw", i, ".Rdata"))
  load(paste0(output_folder, "/output_1_results_draw", i, ".Rdata"))
  load(paste0(output_folder, "/output_3_results_draw", i, ".Rdata"))
  load(paste0(output_folder, "/output_5_results_draw", i, ".Rdata"))
  results_list_base[[i]] <- resbase
  results_list_1[[i]] <- res1
  results_list_3[[i]] <- res3
  results_list_5[[i]] <- res5
}

# RatioDat
SubDat_List<- list()
for (i in 5:7){
  SubDat_List[[i]] <- fread(paste0(project_folder, "/data/orig/SubDataFiles/SubData", Shortyear.Vec[i], "7.csv"))
  print(paste0("Before the procedure, there are ", nrow(SubDat_List[[i]]), "observations"))
}
SubDat <- do.call("rbind", SubDat_List)
SubDat[, group_size := .N, by = c("contractnum", "year")]
RatioDat <- SubDat[, .(
  RatioSum,
  contractnum, 
  famtype_score = fcase(
    nspouse == 0 & ndeps == 0, 1,
    nspouse == 0 & ndeps > 0, 1.85,
    nspouse > 0 & ndeps == 0, 2,
    nspouse > 0 & ndeps > 0, 2.85),
  subscriberid, 
  year), 
]

# Load household - option datasets
ravec = 1:7
yrvec = 2015:2016
nu = .65
source(paste0(project_folder, "/library/DA01QuanModifyExplodedDFunc.R"))
explodeddata_sg_pre <- fread(paste0(project_folder, "/data/explodeddata_sg_pre.csv"))
explodeddata_sg_pre <- explodeddata_sg_pre[ year %in% yrvec]

# Creating final dataset
acgquinsvec <- as.numeric(scan(file=paste0(project_folder, "/data/orig/acg_positive_quintiles"), sep = ",", what = "character")[2:7])
mod_exploded_data <- mod_exploded_data(
  explodeddata_sg_pre, 
  paste0(project_folder, "/data"),
  highmeanthresh = 0.59616,
  lowsumthresh = 0.2575,
  highsumthresh = 1.76329,
  highmaxthresh = 3.37157,
  medincome = 2.746942,
  acgquins = acgquinsvec,
  winsorbound = 89.3707340916662,
  ravec = ravec,
  yrvec = yrvec)

# Kondo.
rm(explodeddata_sg_pre)

# Adjusting to reflect SG institutional details
mod_exploded_data$silver_sub <- 0
mod_exploded_data$silver_costshare <- 0
mod_exploded_data$x_av <- ifelse(
  mod_exploded_data$x_av %in% c(0.7, 0.73, 0.87, 0.94), 
  0.7, 
  mod_exploded_data$x_av)
  
# Adding variables used to create "tiered prices"
mod_exploded_data <- merge(
  mod_exploded_data, 
  RatioDat, 
  by = c("subscriberid", "year"), 
  all.x = T)

results_list_list <- list(
  results_list_1,
  NA,
  results_list_3,
  NA,
  results_list_5)

# Create table of end parameters
make_table <- function(results_list_list0, res_nmb){
  table_list <- list()
  for(i in seq_along(results_list_list0[[res_nmb]])){
    print(paste0("Bootstrap draw ", i, " table is generating.")) 
    table_list[[i]] <- data.table(
      im.post.cs_1 = results_list_list0[[res_nmb]][[i]][[1]]$welfare["CS_1"],
      im.post.cs_2 = results_list_list0[[res_nmb]][[i]][[1]]$welfare["CS_2"],
      im.post.g_spend = results_list_list0[[res_nmb]][[i]][[1]]$welfare["GovmntSpend"], 
      im.post.unin = results_list_list0[[res_nmb]][[i]][[1]]$mktshr["PUninsured"], 
      im.post.bron = results_list_list0[[res_nmb]][[i]][[1]]$mktshr["PBronze"], 
      im.post.silv = results_list_list0[[res_nmb]][[i]][[1]]$mktshr["PSilver"], 
      im.post.gold = results_list_list0[[res_nmb]][[i]][[1]]$mktshr["PGold"], 
      im.post.gross_p_b = results_list_list0[[res_nmb]][[i]][[1]]$grossprem[1], 
      im.post.gross_p_s = results_list_list0[[res_nmb]][[i]][[1]]$grossprem[2],
      im.post.gross_p_g = results_list_list0[[res_nmb]][[i]][[1]]$grossprem[3],
      im.post.net_p_b = results_list_list0[[res_nmb]][[i]][[1]]$netprem[1],
      im.post.net_p_s = results_list_list0[[res_nmb]][[i]][[1]]$netprem[2],
      im.post.net_p_g = results_list_list0[[res_nmb]][[i]][[1]]$netprem[3],
      im.pre.cs_1 = results_list_base[[i]][[1]]$welfare["CS_1"],
      im.pre.cs_2 = results_list_base[[i]][[1]]$welfare["CS_2"],
      im.pre.g_spend = results_list_base[[i]][[1]]$welfare["GovmntSpend"], 
      im.pre.unin = results_list_base[[i]][[1]]$mktshr["PUninsured"], 
      im.pre.bron = results_list_base[[i]][[1]]$mktshr["PBronze"], 
      im.pre.silv = results_list_base[[i]][[1]]$mktshr["PSilver"], 
      im.pre.gold = results_list_base[[i]][[1]]$mktshr["PGold"], 
      im.pre.gross_p_b = results_list_base[[i]][[1]]$grossprem[1], 
      im.pre.gross_p_s = results_list_base[[i]][[1]]$grossprem[2],
      im.pre.gross_p_g = results_list_base[[i]][[1]]$grossprem[3],
      im.pre.net_p_b = results_list_base[[i]][[1]]$netprem[1],
      im.pre.net_p_s = results_list_base[[i]][[1]]$netprem[2],
      im.pre.net_p_g = results_list_base[[i]][[1]]$netprem[3],
      sg.post.cs_1 = results_list_list0[[res_nmb]][[i]][[2]]$welfare["CS_1"],
      sg.post.cs_2 = results_list_list0[[res_nmb]][[i]][[2]]$welfare["CS_2"],
      sg.post.g_spend = results_list_list0[[res_nmb]][[i]][[2]]$welfare["GovmntSpend"], 
      sg.post.unin = results_list_list0[[res_nmb]][[i]][[2]]$mktshr["PUninsured"], 
      sg.post.bron = results_list_list0[[res_nmb]][[i]][[2]]$mktshr["PBronze"], 
      sg.post.silv = results_list_list0[[res_nmb]][[i]][[2]]$mktshr["PSilver"], 
      sg.post.gold = results_list_list0[[res_nmb]][[i]][[2]]$mktshr["PGold"], 
      sg.post.gross_p_b = results_list_list0[[res_nmb]][[i]][[2]]$grossprem[1], 
      sg.post.gross_p_s = results_list_list0[[res_nmb]][[i]][[2]]$grossprem[2],
      sg.post.gross_p_g = results_list_list0[[res_nmb]][[i]][[2]]$grossprem[3],
      sg.post.net_p_b = results_list_list0[[res_nmb]][[i]][[2]]$netprem[1],
      sg.post.net_p_s = results_list_list0[[res_nmb]][[i]][[2]]$netprem[2],
      sg.post.net_p_g = results_list_list0[[res_nmb]][[i]][[2]]$netprem[3],
      sg.pre.cs_1 = results_list_base[[i]][[2]]$welfare["CS_1"],
      sg.pre.cs_2 = results_list_base[[i]][[2]]$welfare["CS_2"],
      sg.pre.g_spend = results_list_base[[i]][[2]]$welfare["GovmntSpend"], 
      sg.pre.unin = results_list_base[[i]][[2]]$mktshr["PUninsured"], 
      sg.pre.bron = results_list_base[[i]][[2]]$mktshr["PBronze"], 
      sg.pre.silv = results_list_base[[i]][[2]]$mktshr["PSilver"], 
      sg.pre.gold = results_list_base[[i]][[2]]$mktshr["PGold"], 
      sg.pre.gross_p_b = results_list_base[[i]][[2]]$grossprem[1], 
      sg.pre.gross_p_s = results_list_base[[i]][[2]]$grossprem[2],
      sg.pre.gross_p_g = results_list_base[[i]][[2]]$grossprem[3],
      sg.pre.net_p_b = results_list_base[[i]][[2]]$netprem[1],
      sg.pre.net_p_s = results_list_base[[i]][[2]]$netprem[2],
      sg.pre.net_p_g = results_list_base[[i]][[2]]$netprem[3]
    )
  }
  table <- do.call("rbind", table_list)
  
  print("generating small group pre surplus:")
  for(i in seq_along(results_list_list0[[res_nmb]])){
    print(paste0("Bootstrap draw ", i, " SG pre surplus is being estimated.")) 
    # Adding in small group consumer surplus pre
    ## Load small group parameters
    SGMarketEstimates <- fread(paste0(project_folder, "/analysis/estimationcode/specs/lowrisk_simplemh_censor0.2_switcher/bootstrap/output_", i, "/solution.csv"))
    SGMarketEstimates <- SGMarketEstimates[, -c(1,2)]
    SGMarketlabels  <- fread(paste0(project_folder, "/analysis/estimationcode/specs/lowrisk_simplemh_censor0.2_switcher/bootstrap/output_", i, "/covar_lens.csv"))
    SGMarSE <- cumsum(SGMarketlabels$x)
    ## Create the small group surplus calculation
    source(paste0(project_folder, "/analysis/counterfactuals/bootstrap/FA01CreateSmallGroupSurplusDataset.R"), local = T)
    ## dropping outside options and keeping 2016
    sg_surplus <- sg_surplus[!(x_av == 0)]
    sg_surplus <- sg_surplus[year == 2016]
    sg_surplus <- sg_surplus[!(fixedeffect_beta0 == 0 & x_av != 0)]
    # Generating surplus:
    sg_sur_collapsed <- sg_surplus[
      x_av > 0, 
      .(logsum = log(sum(exp(V))), alphampsi = mode1(alpha - psi)), 
      by = c("subscriberid", "year")
    ]
    sg_sur_collapsed$cs <- sg_sur_collapsed$logsum / sg_sur_collapsed$alphampsi
    table$sg.pre.cs_1[i] <- mean(sg_sur_collapsed[alphampsi > alphampsicutoff]$cs)
    table$sg.pre.cs_2[i] <- mean(sg_sur_collapsed[alphampsi > alphampsicutoff]$cs)
    print(paste0(i))
  }
  
  # Create change in small group CS:
  table$sg.change.cs_1 <- table$sg.post.cs_1 - table$sg.pre.cs_1
  table$sg.change.cs_2 <- table$sg.post.cs_2 - table$sg.pre.cs_2
  table$im.change.cs_1 <- table$im.post.cs_1 - table$im.pre.cs_1
  table$im.change.cs_2 <- table$im.post.cs_2 - table$im.pre.cs_2
  return(table)
}

# save as
table_res1 <- make_table(results_list_list, 1)
fwrite(table_res1, paste0(output_folder, "/outcome_table_res1.csv"))
table_res3 <- make_table(results_list_list, 3)
fwrite(table_res3, paste0(output_folder, "/outcome_table_res3.csv"))
table_res5 <- make_table(results_list_list, 5)
fwrite(table_res5, paste0(output_folder, "/outcome_table_res5.csv"))

# create the final tables 
column_vector_base <- c("im.change.cs_2", "im.pre.g_spend", "im.pre.unin", "im.pre.bron", "im.pre.silv", "im.pre.gold",  
  "im.pre.gross_p_b", "im.pre.gross_p_s", "im.pre.gross_p_g", "im.pre.net_p_b", "im.pre.net_p_s", "im.pre.net_p_g", 
  "sg.change.cs_2", "sg.pre.g_spend", "sg.pre.unin", "sg.pre.bron", "sg.pre.silv", "sg.pre.gold",  
  "sg.pre.gross_p_b", "sg.pre.gross_p_s", "sg.pre.gross_p_g", "sg.pre.net_p_b", "sg.pre.net_p_s", "sg.pre.net_p_g")
  column_vector_post <- c("im.change.cs_2", "im.post.g_spend", "im.post.unin", "im.post.bron", "im.post.silv", "im.post.gold",  
  "im.post.gross_p_b", "im.post.gross_p_s", "im.post.gross_p_g", "im.post.net_p_b", "im.post.net_p_s", "im.post.net_p_g", 
  "sg.change.cs_2", "sg.post.g_spend", "sg.post.unin", "sg.post.bron", "sg.post.silv", "sg.post.gold",  
  "sg.post.gross_p_b", "sg.post.gross_p_s", "sg.post.gross_p_g", "sg.post.net_p_b", "sg.post.net_p_s", "sg.post.net_p_g")
base_vec <- lapply(table_res1[, column_vector_base, with = F], sd)
res1_vec <- lapply(table_res1[, column_vector_post, with = F], sd)
res3_vec <- lapply(table_res3[, column_vector_post, with = F], sd)
res5_vec <- lapply(table_res5[, column_vector_post, with = F], sd)

results_mat <- data.table(
  names = names(unlist(base_vec)), 
  base_vec, 
  res1_vec, 
  res3_vec, 
  res5_vec)
  
fwrite(results_mat, paste0(project_folder, "/analysis/tablesandfigures/release/counterfactuals/bootstrap/bootstrap_stderr.csv"))
